##############################
####
#### Authors: Stefan Müller, and Tom Louwerse
#### Title: "The Electoral Cycle Effect in Parliamentary Democracies."
#### Journal: Political Science Research and Methods
#### 
#### Replication Material

##############################


### File description: This file reproduces all regressions, plots, 
### and tables reported in the paper and the Supporting Information. 
### Note that the NLME models can take a long time to run. 
### Therefore, we moved these models, tables, and plots 
### into a separate R script (02_models_plots_tables_nlme.R).


# 1) Load libraries and prepare data -----

library(tidyverse)
library(lme4)
library(effects)
library(texreg)
library(car)
library(pastecs)
library(stargazer)
library(ggthemes)
library(nlme)

# Load functions for texreg tables and ggplot2 theme
source("helper.R")

# Set theme for all plots
theme_set(theme_base_edited())

# Load data
dta <- readRDS("data.rds")

# Apply filters
dta_all <- dta %>%
  dplyr::filter(total_length >= 182.2) %>%
  dplyr::filter(electoral_cycle_cabinet > 0, electoral_cycle_max < 1.2) %>%
  dplyr::filter(election == "Legislative") %>%
  dplyr::filter(system != "Presidential") %>% 
  mutate(country = as.factor(country),
         country_code = as.factor(country_code),
         election_id = as.factor(election_id),
         poll_id = as.factor(poll_id)) %>% 
    mutate(partyid = as.factor(paste(country, partyid, sep = "_")))

# Only select polls by government parties
dta_gov <- dta_all %>%
  dplyr::filter(government == 1 & electoral_cycle_cabinet > 0)


# Removing missing values required to plot resid vs fitted
dta_gov_base <- dta_gov %>% 
  filter(!is.na(poll_change) & !is.na(electoral_cycle_cabinet)
         & !is.na(gdp_change_lag) & !is.na(vote_share_last) & !is.na(elecyear))

### 2) Run regression models ----

## Basic model
reg_base <- lmer(poll_change ~ 
                   poly(electoral_cycle_cabinet, 3)  + 
                   gdp_change_lag +
                   vote_share_last + 
                   I(elecyear - 1986) + 
                   (1 | country) +  
                   (1 |  country: partyid) + 
                   (1 | election_id : partyid), data = dta_gov_base)


# Removing missing values (required to plot resid vs fitted)
dta_gov_base <- dta_gov %>% 
  filter(!is.na(poll_change) & !is.na(electoral_cycle_cabinet)
         & !is.na(gdp_change_lag) & !is.na(vote_share_last) & !is.na(elecyear)) %>% 
  group_by(country, election_id, partyid) %>%
  arrange(country, election_id, partyid, polldate)


## Basic model without cubic term
reg_base_square <- lmer(poll_change ~ 
                          poly(electoral_cycle_cabinet, 2)  + 
                          gdp_change_lag +
                          vote_share_last + 
                          I(elecyear - 1986) + 
                          (1 | country) +  
                          (1 | country: partyid) + 
                          (1 | election_id : partyid), data = dta_gov_base)


# Removing missing values (required to plot resid vs fitted)
dta_gov_base_lag <- dta_gov_base %>% 
  filter(!is.na(poll_change_lag))

## Basic model with lag
reg_base_lag <- update(reg_base, . ~ . + poll_change_lag, data = dta_gov_base_lag) 


## Combined model

# Removing missing values (required to plot resid vs fitted)
dta_gov_base_combined <- dta_gov_base %>% 
  filter(!is.na(diss_pm) & !is.na(single_party_cabinet))

reg_combined <- update(reg_base, . ~ . + poly(electoral_cycle_cabinet, 3) * diss_pm +
                         poly(electoral_cycle_cabinet, 3) * single_party_cabinet, 
                       data = dta_gov_base_combined) 

## Combined model (and only countries with three or more cycles)

countries_remove <- c("Belgium", "Finland", "Hungary", 
                      "Iceland", "Italy", "Poland",
                      "Switzerland", "Turkey")

dta_gov_at_leat_three_cycles <- dta_gov %>% 
  filter(!country_code %in% countries_remove)

reg_combined_subset <- update(reg_base, . ~ . + poly(electoral_cycle_cabinet, 3) * diss_pm +
                                poly(electoral_cycle_cabinet, 3) * single_party_cabinet,
                              data = dta_gov_at_leat_three_cycles) # and another 

## Use percentage change as DV
reg_base_percentage <- lmer(poll_change_percentage_av ~ 
                              poly(electoral_cycle_cabinet, 3)  + 
                              gdp_change_lag +
                              vote_share_last + 
                              I(elecyear - 1986) + 
                              (1|country) +  
                              (1 |  country: partyid) + 
                              (1 | election_id : partyid), data = dta_gov) 


## Single-party cabinet
reg_single_party_gov <- update(reg_base, . ~ . + 
                                 poly(electoral_cycle_cabinet, 3) * single_party_cabinet) 


## Prime minister dissolution power
reg_pm_diss <- update(reg_base, . ~ . + poly(electoral_cycle_cabinet, 3) * diss_pm) 


## Government dissolution power
reg_gov_diss <- update(reg_base, . ~ . + poly(electoral_cycle_cabinet, 3) * diss_gov) 

## Minority cabinet
reg_min_gov <- update(reg_base, . ~ . + poly(electoral_cycle_cabinet, 3) * minority_cabinet)

## Largest government party
dta_gov_coa <- subset(dta_gov, single_party_cabinet==0)
reg_largest_gov_party <- update(reg_base, . ~ . + 
                                  poly(electoral_cycle_cabinet, 3) * largest_government_party,
                                data = dta_gov_coa)
## Lijphart scores
reg_lijphart <- update(reg_base, . ~ . + poly(electoral_cycle_cabinet, 3) * exec_parties_4510 +
                         poly(electoral_cycle_cabinet, 3) *  fed_unit_4510) 


## Decades

# Select countries where we have polling data since the 1960s
country_selection <- c("Australia", "Canada", 
                       "Denmark", "Germany", 
                       "Netherlands", "Norway", 
                       "United Kingdom")

dta_decades <- subset(dta_gov, elecyear > 1959 & 
                        country_code %in% country_selection)

dta_decades <- dta_decades %>% 
  dplyr::arrange(country, elecdecade) %>% 
  dplyr::filter(elecyear > "1959")

model_decades <- lmer(poll_change ~ poly(electoral_cycle_cabinet, 3) * elecdecade + 
                        gdp_change_lag +
                        vote_share_last + 
                        (1| country) +  
                        (1 | country: partyid) + 
                        (1 | election_id : partyid), data = dta_decades)


reg_combined_gdp <- update(reg_combined, . ~  
                             poly(electoral_cycle_cabinet, 3) * gdp_change_lag + ., 
                           data = dta_gov_base_combined)


## Standardized GDP change (for each country and decade)
reg_combined_gdp_stand <- update(reg_combined, . ~  
                                   poly(electoral_cycle_cabinet, 3) * gdp_change_lag_stand + ., 
                                 data = dta_gov_base_combined)

## Alternative measurement of electoral_cycle
reg_combined_actual_cycle <- lmer(poll_change ~ 
                                    poly(electoral_cycle_planned, 3)  +
                                    poly(electoral_cycle_planned, 3) * diss_pm + 
                                    poly(electoral_cycle_planned, 3) * single_party_cabinet +
                                    gdp_change_lag +
                                    vote_share_last + I(elecyear - 1986) + 
                                    (1|country) + 
                                    (1 |  country: partyid) + 
                                    (1 | election_id : partyid), data = filter(dta_gov, electoral_cycle_planned <= 1))


## Type of election
reg_combined_elec_type <- update(reg_combined, . ~  
                                   poly(electoral_cycle_cabinet, 3) * early_election_type +
                                   single_party_cabinet + ., 
                                 data = dta_gov)


#### 3) Create plots from paper -----

## Fig 2a: Effect plot: Basic model ----
eff_base <- data.frame(Effect("electoral_cycle_cabinet", reg_base,
                              xlevels = 30))

ggplot(eff_base, aes(x = electoral_cycle_cabinet,
                     y = fit, ymin = lower, ymax = upper)) + 
    geom_line() +
    geom_ribbon(fill = "grey50", alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change")


## Fig. 2b: Effect plot: Single party government -----
eff_single_party_gov <- data.frame(Effect(c("electoral_cycle_cabinet", "single_party_cabinet"), 
                                          reg_single_party_gov,
                                          xlevels= 30))

ggplot(eff_single_party_gov, aes(x = electoral_cycle_cabinet,
                                 y = fit, ymin = lower, ymax = upper, 
                                 fill = factor(single_party_cabinet), 
                                 linetype = factor(single_party_cabinet))) + 
    geom_line() +
    geom_ribbon(alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    scale_fill_manual(values = c("grey60", "dodgerblue2"), 
                      name = "Single Party Government", 
                      labels = c("No", "Yes"),
                      guide = guide_legend(reverse = TRUE)) +
    scale_linetype_manual(values = c(1, 2), 
                          name = "Single Party Government", 
                          labels = c("No", "Yes"),
                          guide = guide_legend(reverse = TRUE)) +
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change") +
    theme(legend.justification = c(0, 0),
          legend.position = c(0.4, 0.7))

## Fig 2c: Effect plot: PM dissolution power ----
eff_pm_diss <- data.frame(Effect(c("electoral_cycle_cabinet", "diss_pm"), reg_pm_diss,
                                 xlevels=list(diss_pm=c(0,10), 
                                              electoral_cycle_cabinet = seq(0,1,0.01))))

ggplot(eff_pm_diss, aes(x = electoral_cycle_cabinet,
                        y = fit, ymin = lower, ymax = upper, 
                        fill = factor(diss_pm), 
                        linetype = factor(diss_pm))) + 
    geom_line() +
    geom_ribbon(alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    scale_fill_manual(values = c("grey60", "dodgerblue2"), 
                      name = "PM Dissolution Power", 
                      labels = c("Lowest (0)", "Highest (10)"),
                      guide = guide_legend(reverse = TRUE)) +
    scale_linetype_manual(values = c(1, 2), 
                          name = "PM Dissolution Power", 
                          labels = c("Lowest (0)", "Highest (10)"),
                          guide = guide_legend(reverse = TRUE)) +
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change") +
    theme(legend.justification = c(0, 0),
          legend.position = c(0.55, 0.7))



## Fig. 3: Historical development of the electoral cycle effect ----

effect_decade <- as.data.frame(Effect(c("electoral_cycle_cabinet", "elecdecade"),
                                      model_decades,
                                      xlevels=list(elecdecade=unique(dta_decades$elecdecade), 
                                                   electoral_cycle_cabinet=seq(0,1,.01)))) %>% 
    mutate(elecdecade_print = paste("Decade:", elecdecade))

ggplot(data = data.frame(effect_decade), aes(x = electoral_cycle_cabinet, y = fit)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey80") +
    geom_line() + 
    facet_wrap(~elecdecade_print) +
    scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
    xlab("Electoral Cycle") +
    ylab("Expected Absolute Poll Change") 


##### 4) Create plots from appendix ----

## Fig. A1: Loess regression line with 95 per cent confidence intervals
## (takes some time to run) ----
dta_gov_cab_cycle <- subset(dta_gov, electoral_cycle_cabinet > 0)

ggplot(data = dta_gov_cab_cycle, aes(x = electoral_cycle_cabinet, y = poll_change)) +
    geom_smooth(method = "loess", se = TRUE, colour = "black") +
    geom_hline(yintercept = 0, linetype = 2, colour = "black") +
    xlab("Electoral Cycle") +
    ylab("Change in Party Support")


## Note: Fig. A2a, A2b, A2c, and A3 are part of the file models_plots_tables_nlme.R

## Fig. A4a: Basic model ----
ggplot(reg_base, aes(.fitted, .resid)) + 
    geom_point(alpha = 0.3) +
    geom_smooth(n = 40) +
    labs(x = "Fitted values", y = "Residuals")

## Fig. A4b: Combined Model -----   
ggplot(reg_combined, aes(.fitted, .resid)) + 
    geom_point(alpha = 0.3) +
    geom_smooth() +
    labs(x = "Fitted values", y = "Residuals")

## Fig. A4c: Basic model with lagged DV ----
ggplot(reg_base_lag, aes(.fitted, .resid)) + 
    geom_point(alpha = 0.3) +
    geom_smooth() +
    labs(x = "Fitted values", y = "Residuals")


# Fig A5a: Effect plot: Minority government ----
eff_min_gov <- data.frame(Effect(c("electoral_cycle_cabinet", "minority_cabinet"), 
                                 reg_min_gov,
                                 xlevels= 30))

ggplot(eff_min_gov, aes(x = electoral_cycle_cabinet,
                        y = fit, ymin = lower, ymax = upper, 
                        fill = factor(minority_cabinet), 
                        linetype = factor(minority_cabinet))) + 
    geom_line() +
    geom_ribbon(alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    scale_fill_manual(values = c("grey60", "dodgerblue2"), 
                      name = "Minority Government", 
                      labels = c("No", "Yes")) +
    scale_linetype_manual(values = c(1, 2), 
                          name = "Minority Government", 
                          labels = c("No", "Yes")) +
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change") +
    theme(legend.justification = c(0, 0),
          legend.position = c(0.55, 0.7))

## Fig. A5b: Effect plot: Largest coalition government ----
eff_largest_gov <- data.frame(Effect(c("electoral_cycle_cabinet", "largest_government_party"), 
                                     reg_largest_gov_party,
                                     xlevels= 30))

ggplot(eff_largest_gov, aes(x = electoral_cycle_cabinet,
                            y = fit, ymin = lower, ymax = upper, 
                            fill = factor(largest_government_party), 
                            linetype = factor(largest_government_party))) + 
    geom_line() +
    geom_ribbon(alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    scale_fill_manual(values = c("grey60", "dodgerblue2"), name = "Largest Government Party", 
                      labels = c("No", "Yes"),
                      guide = guide_legend(reverse = TRUE)) +
    scale_linetype_manual(values = c(1, 2), 
                          name = "Largest Government Party", 
                          labels = c("No", "Yes"),
                          guide = guide_legend(reverse = TRUE)) +
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change") +
    theme(legend.justification = c(0, 0),
          legend.position = c(0.4, 0.7))


## Fig. A5c: Executive parties dimension ----
dta_effect_exec_parties <- as.data.frame(Effect(c("electoral_cycle_cabinet", "exec_parties_4510"),
                                                reg_lijphart, 
                                                xlevels=list(exec_parties_4510 = c(as.numeric(min(dta_gov$exec_parties_4510, na.rm = TRUE)), 
                                                                                   as.numeric(mean(dta_gov$exec_parties_4510, na.rm = TRUE)),
                                                                                   as.numeric(max(dta_gov$exec_parties_4510, na.rm = TRUE))), 
                                                             electoral_cycle_cabinet=seq(0,1,.01)))) %>% 
    mutate(exec_parties_print = ifelse(exec_parties_4510 == as.numeric(min(dta_gov$exec_parties_4510, na.rm = TRUE)), 'Executive: Strong', 
                                       ifelse(exec_parties_4510 > as.numeric(min(dta_gov$exec_parties_4510, na.rm = TRUE)) & exec_parties_4510 < as.numeric(max(dta_gov$exec_parties_4510, na.rm = TRUE)), 'Medium', 
                                              'Executive: Weak')))

dta_effect_exec_parties$exec_parties_print <- factor(dta_effect_exec_parties$exec_parties_print, 
                                                     levels = c("Executive: Strong", "Medium", "Executive: Weak"))

ggplot(data = dta_effect_exec_parties, aes(x = electoral_cycle_cabinet, y = fit)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey80") +
    geom_line() + 
    facet_wrap(~exec_parties_print) +
    scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
    xlab("Electoral Cycle") +
    ylab("Expected Absolute Poll Change") 

## Fig A5d: Federal unitary dimension ----
dta_effect_fed_unitary <- as.data.frame(Effect(c("electoral_cycle_cabinet", "fed_unit_4510"),
                                               reg_lijphart, 
                                               xlevels=list(fed_unit_4510=c(as.numeric(min(dta_gov$fed_unit_4510, na.rm = TRUE)), 
                                                                            as.numeric(mean(dta_gov$fed_unit_4510, na.rm = TRUE)),
                                                                            as.numeric(max(dta_gov$fed_unit_4510, na.rm = TRUE))), 
                                                            electoral_cycle_cabinet=seq(0,1,.01)))) %>% 
    mutate(fed_unit_print = ifelse(fed_unit_4510 == as.numeric(min(dta_gov$fed_unit_4510, na.rm = TRUE)), 'Unitary', 
                                   ifelse(fed_unit_4510 > as.numeric(min(dta_gov$fed_unit_4510, na.rm = TRUE)) & fed_unit_4510 < as.numeric(max(dta_gov$fed_unit_4510, na.rm = TRUE)), 'Medium', 
                                          'Federal')))

dta_effect_fed_unitary$fed_unit_print <- factor(dta_effect_fed_unitary$fed_unit_print, 
                                                levels = c("Unitary", "Medium", "Federal"))

ggplot(data = dta_effect_fed_unitary, aes(x = electoral_cycle_cabinet, y = fit)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey80") +
    geom_line() + 
    facet_wrap(~fed_unit_print) +
    scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
    xlab("Electoral Cycle") +
    ylab("Expected Absolute Poll Change") 


## Fig. A5e: Effect plot: Government dissolution power -----
eff_gov_diss <- data.frame(Effect(c("electoral_cycle_cabinet", "diss_gov"), reg_gov_diss,
                                  xlevels=list(diss_gov=c(0,8.5), 
                                               electoral_cycle_cabinet = seq(0,1,0.01))))

ggplot(eff_gov_diss, aes(x = electoral_cycle_cabinet,
                         y = fit, ymin = lower, ymax = upper, 
                         fill = factor(diss_gov), 
                         linetype = factor(diss_gov))) + 
    geom_line() +
    geom_ribbon(alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    scale_fill_manual(values = c("grey60", "dodgerblue2"), name = "Government Dissolution Power", 
                      labels = c("Lowest (0)", "Highest (8.5)"),
                      guide = guide_legend(reverse = TRUE)) +
    scale_linetype_manual(values = c(1, 2), 
                          name = "Government Dissolution Power", 
                          labels = c("Lowest (0)", "Highest (8.5)"),
                          guide = guide_legend(reverse = TRUE)) +
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change") +
    theme(legend.justification = c(0, 0),
          legend.position = c(0.4, 0.7))

## Fig. A6a = Fig. 2b ----


## Fig. A6b: Effect plot: Single party government (with combined model) ----
eff_single_party_gov_combined <- data.frame(Effect(c("electoral_cycle_cabinet", "single_party_cabinet"), 
                                                   reg_combined,
                                                   xlevels= 30))

ggplot(eff_single_party_gov_combined, aes(x = electoral_cycle_cabinet,
                                          y = fit, ymin = lower, ymax = upper, 
                                          fill = factor(single_party_cabinet), 
                                          linetype = factor(single_party_cabinet))) + 
    geom_line() +
    geom_ribbon(alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    scale_fill_manual(values = c("grey60", "dodgerblue2"), name = "Single Party Government", 
                      labels = c("No", "Yes"),
                      guide = guide_legend(reverse = TRUE)) +
    scale_linetype_manual(values = c(1, 2), 
                          name = "Single Party Government", 
                          labels = c("No", "Yes"),
                          guide = guide_legend(reverse = TRUE)) +
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change") +
    theme(legend.justification = c(0, 0),
          legend.position = c(0.4, 0.7))


## Fig. A5c = Fig. 2c ----


## Fig. A6d: Effect plot: PM dissolution power (with combined model) ----
eff_pm_diss_combined <- data.frame(Effect(c("electoral_cycle_cabinet", "diss_pm"), reg_combined,
                                          xlevels=list(diss_pm=c(0,10), 
                                                       electoral_cycle_cabinet = seq(0,1,0.01))))

ggplot(eff_pm_diss_combined, aes(x = electoral_cycle_cabinet,
                                 y = fit, ymin = lower, ymax = upper, 
                                 fill = factor(diss_pm), 
                                 linetype = factor(diss_pm))) + 
    geom_line() +
    geom_ribbon(alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    scale_fill_manual(values = c("grey60", "dodgerblue2"), name = "PM Dissolution Power", 
                      labels = c("Lowest (0)", "Highest (10)"),
                      guide = guide_legend(reverse = TRUE)) +
    scale_linetype_manual(values = c(1, 2), 
                          name = "PM Dissolution Power", 
                          labels = c("Lowest (0)", "Highest (10)"),
                          guide = guide_legend(reverse = TRUE)) +
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change") +
    theme(legend.justification = c(0, 0),
          legend.position = c(0.52, 0.7))


## Fig. A7a = Fig. A2b ----


## Fig. A7b: Effect plot: Single party government (with combined model and more than three cycles per country) ----
eff_single_party_gov_combined_subset <- data.frame(Effect(c("electoral_cycle_cabinet", "single_party_cabinet"), 
                                                          reg_combined_subset,
                                                          xlevels= 30))

ggplot(eff_single_party_gov_combined_subset, aes(x = electoral_cycle_cabinet,
                                                 y = fit, ymin = lower, ymax = upper, 
                                                 fill = factor(single_party_cabinet), 
                                                 linetype = factor(single_party_cabinet))) + 
    geom_line() +
    geom_ribbon(alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    scale_fill_manual(values = c("grey60", "dodgerblue2"), name = "Single Party Government", 
                      labels = c("No", "Yes"),
                      guide = guide_legend(reverse = TRUE)) +
    scale_linetype_manual(values = c(1, 2), 
                          name = "Single Party Government", 
                          labels = c("No", "Yes"),
                          guide = guide_legend(reverse = TRUE)) +
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change") +
    theme(legend.justification = c(0, 0),
          legend.position = c(0.4, 0.7))

## Fig. A7c = Fig. A2c ----


## Fig. A7d: Effect plot: PM dissolution power (with combined model and more than three cycles per country) ----
eff_pm_diss_combined_subset <- data.frame(Effect(c("electoral_cycle_cabinet", "diss_pm"), reg_combined_subset,
                                                 xlevels=list(diss_pm=c(0,10), 
                                                              electoral_cycle_cabinet = seq(0,1,0.01))))

ggplot(eff_pm_diss_combined_subset, aes(x = electoral_cycle_cabinet,
                                        y = fit, ymin = lower, ymax = upper, 
                                        fill = factor(diss_pm), 
                                        linetype = factor(diss_pm))) + 
    geom_line() +
    geom_ribbon(alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    scale_fill_manual(values = c("grey60", "dodgerblue2"), name = "PM Dissolution Power", 
                      labels = c("Lowest (0)", "Highest (10)"),
                      guide = guide_legend(reverse = TRUE)) +
    scale_linetype_manual(values = c(1, 2), 
                          name = "PM Dissolution Power", 
                          labels = c("Lowest (0)", "Highest (10)"),
                          guide = guide_legend(reverse = TRUE)) +
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change") +
    theme(legend.justification = c(0, 0),
          legend.position = c(0.52, 0.7))


## Fig A8a Effect plot: Basic model with lagged DV ----
eff_base_lag <- data.frame(Effect("electoral_cycle_cabinet", reg_base_lag,
                                  xlevels = 30))

ggplot(eff_base_lag, aes(x = electoral_cycle_cabinet,
                         y = fit, ymin = lower, ymax = upper)) + 
    geom_line() +
    geom_ribbon(fill = "grey50", alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change")


## Fig A8b = Fig 2a ----

## Fig A8c: Effect plot: Basic model percentage ----
eff_base_perc <- data.frame(Effect("electoral_cycle_cabinet", reg_base_percentage,
                                   xlevels = 30))

ggplot(eff_base_perc, aes(x = electoral_cycle_cabinet,
                          y = fit, ymin = lower, ymax = upper)) + 
    geom_line() +
    geom_ribbon(fill = "grey50", alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    labs(x = "Electoral Cycle", y = "Expected Poll Change (in %)")


## Fig. A9: GDP change ----
dta_effect_gdp <- as.data.frame(Effect(c("electoral_cycle_cabinet", "gdp_change_lag"),
                                       reg_combined_gdp, 
                                       xlevels=list(gdp_change_lag=c(as.numeric(min(dta_gov_base_combined$gdp_change_lag, na.rm = TRUE)), 
                                                                     as.numeric(mean(dta_gov_base_combined$gdp_change_lag, na.rm = TRUE)),
                                                                     as.numeric(max(dta_gov_base_combined$gdp_change_lag, na.rm = TRUE))), 
                                                    electoral_cycle_cabinet=seq(0,1,.01)))) %>% 
    mutate(gdp_change_lag_print = ifelse(gdp_change_lag == as.numeric(min(dta_gov_base_combined$gdp_change_lag, na.rm = TRUE)), 'GDP Change: Low', 
                                         ifelse(gdp_change_lag > as.numeric(min(dta_gov_base_combined$gdp_change_lag, na.rm = TRUE)) & gdp_change_lag < as.numeric(max(dta_gov_base_combined$gdp_change_lag, na.rm = TRUE)), 'GDP Change: Medium', 
                                                'GDP Change: High')))

dta_effect_gdp$gdp_change_lag_print <- factor(dta_effect_gdp$gdp_change_lag_print, 
                                              levels = c("GDP Change: Low", "GDP Change: Medium", "GDP Change: High"))

ggplot(data = dta_effect_gdp, aes(x = electoral_cycle_cabinet, y = fit)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey80") +
    geom_line() + 
    facet_wrap(~gdp_change_lag_print) +
    scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
    xlab("Electoral Cycle") +
    ylab("Expected Absolute Poll Change") 


## Fig. A10: GDP change (standardised) ----
dta_effect_gdp_stand <- as.data.frame(Effect(c("electoral_cycle_cabinet", "gdp_change_lag_stand"),
                                             reg_combined_gdp_stand, 
                                             xlevels=list(gdp_change_lag_stand=c(as.numeric(min(dta_gov_base_combined$gdp_change_lag_stand, na.rm = TRUE)), 
                                                                                 as.numeric(mean(dta_gov_base_combined$gdp_change_lag_stand, na.rm = TRUE)),
                                                                                 as.numeric(max(dta_gov_base_combined$gdp_change_lag_stand, na.rm = TRUE))), 
                                                          electoral_cycle_cabinet=seq(0,1,.01)))) %>% 
    mutate(gdp_change_lag_stand_print = ifelse(gdp_change_lag_stand == as.numeric(min(dta_gov_base_combined$gdp_change_lag_stand, na.rm = TRUE)), 'GDP Change: Low', 
                                               ifelse(gdp_change_lag_stand > as.numeric(min(dta_gov_base_combined$gdp_change_lag_stand, na.rm = TRUE)) & gdp_change_lag_stand < as.numeric(max(dta_gov_base_combined$gdp_change_lag_stand, na.rm = TRUE)), 'GDP Change: Medium', 
                                                      'GDP Change: High')))

dta_effect_gdp_stand$gdp_change_lag_stand_print <- factor(dta_effect_gdp_stand$gdp_change_lag_stand_print, 
                                                          levels = c("GDP Change: Low", "GDP Change: Medium", "GDP Change: High"))

ggplot(data = dta_effect_gdp_stand, aes(x = electoral_cycle_cabinet, y = fit)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey80") +
    geom_line() + 
    facet_wrap(~gdp_change_lag_stand_print) +
    scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
    xlab("Electoral Cycle") +
    ylab("Expected Absolute Poll Change") 


# Fig. A11: Plot with early/regular/failure elections ----
dta_election_type <- dta %>% 
    dplyr::select(early_election_type, country_code, elecyear, election_id) %>% 
    unique() %>% 
    group_by(country_code, early_election_type) %>% 
    dplyr::summarise(number_elections = n())

ggplot(dta_election_type, aes(x = early_election_type, y = number_elections)) +
    geom_bar(stat = "identity", width = 0.5) +
    facet_wrap(~country_code) +
    coord_flip() + 
    geom_text(aes(label = number_elections), hjust = -0.2, size = 3.5) + 
    labs(y = "Number of elections", x = "Election type") +
    scale_y_continuous(limits = c(0, 22), breaks = c(seq(0, 20, 5))) 


## Fig. A12: Election type ----
dta_effect_early_election_type <- as.data.frame(Effect(c("electoral_cycle_cabinet", "early_election_type"),
                                                       reg_combined_elec_type,
                                                       xlevels=list(electoral_cycle_cabinet=seq(0,1,.01)))) %>% 
    mutate(early_election_type_print = paste("Type:", early_election_type, sep = " "))

dta_effect_early_election_type$early_election_type_print <- factor(dta_effect_early_election_type$early_election_type_print, levels = c("Type: Regular", "Type: Opportunistic", "Type: Failure"))

ggplot(data = dta_effect_early_election_type, aes(x = electoral_cycle_cabinet, y = fit)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey80") +
    geom_line() + 
    facet_wrap(~early_election_type_print) +
    scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
    xlab("Electoral Cycle") +
    ylab("Expected Absolute Poll Change") 



## Fig A13a: Effect plot: Basic model with a different electoral cycle operationalisation  ----
eff_base_combined_planned <- data.frame(Effect("electoral_cycle_cabinet", reg_combined,
                                               xlevels = 30))

ggplot(eff_base_combined_planned, aes(x = electoral_cycle_cabinet,
                                      y = fit, ymin = lower, ymax = upper)) + 
    geom_line() +
    geom_ribbon(fill = "grey50", alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    labs(x = "Electoral Cycle (actual length)", y = "Expected Absolute Poll Change")



## Fig A13b: Effect plot: Basic model with a different electoral cycle operationalisation  ----
eff_base_combined_actual_cycle <- data.frame(Effect("electoral_cycle_planned", reg_combined_actual_cycle,
                                                    xlevels = 30))

ggplot(eff_base_combined_actual_cycle, aes(x = electoral_cycle_planned,
                                           y = fit, ymin = lower, ymax = upper)) + 
    geom_line() +
    geom_ribbon(fill = "grey50", alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    labs(x = "Electoral Cycle (planned length)", y = "Expected Absolute Poll Change")

### 5) Reproduce tables ----


# Recode variable names to nice formats
nice.varnames <- "'poll_change_lag'='Poll change$_{t-1}$';
'poly(electoral_cycle_cabinet, 2)1'='El. cycle';
'poly(electoral_cycle_cabinet, 2)2'='El. cycle$^2$';
'poly(electoral_cycle_cabinet, 3)1'='El. cycle';
'poly(electoral_cycle_cabinet, 3)2'='El. cycle$^2$';
'poly(electoral_cycle_cabinet, 3)3'='El. cycle$^3$';
'poly(electoral_cycle_planned, 3)1'='El. cycle (planned)';
'poly(electoral_cycle_planned, 3)2'='El. cycle$^2$ (planned)';
'poly(electoral_cycle_planned, 3)3'='El. cycle$^3$ (planned)';
'gdp_change_lag'='GDP Change';
'gdp_change_lag_stand'='GDP Change (standardized)';
'unemployment'='Unemployment';
'vote_share_last'='Party support at last election';
'single_party_cabinet1'='Single party gov.';
'minority_cabinet1'='Minority gov.';
'diss_pm'='PM diss. power';
'diss_gov'='Gov. diss. power';
'early_election_typeOpportunistic'='El. --  Opportunistic';
'early_election_typeRegular'='El. -- Regular';
'largest_government_party1'='Largest gov. party';
'I(elecyear - 1986)'='Election year - 1986';
'elecdecade1970'='Decade -- 1970s';
'elecdecade1980'='Decade -- 1980s';
'elecdecade1990'='Decade -- 1990s';
'elecdecade2000'='Decade -- 2000s';
'elecdecade2010'='Decade -- 2010s';
'exec_parties_4510'='Executive-parties dim.';
'fed_unit_4510'='Federal-unitary dim.'"


# Code which makes sure that interaction terms
# are also recoded in to nice variable names
splitRecode <- function(x, by=":", recodes) {
    if(length(x)>1) {
        out <- sapply(x, function(q) splitRecode(q, by=by, recodes=recodes))
        names(out) <- NULL
        return(out)
    }
    y <- unlist(strsplit(x, ":"))
    z <- car::recode(y, nice.varnames)
    if(length(z)>1) z <- paste(z, collapse=" $\\times$ ")
    return(z)
}


## Tab A1: Descriptive Statistics

nice.names <-c("Poll change",
               "Poll change$_{t-1}$",
               "Electoral cycle (planned)",
               "Electoral cycle (cabinet)",
               "Single party government",
               "Largest government party",
               "Minority government",
               "GDP change",
               "GDP change (standardized)",
               "Election: Regular",
               "Election: Failure",
               "Election: Opportunistic",
               "PM dissolution power",
               "Government dissolution power",
               "Party support at last election")


desc_table <- dta_gov %>%
    filter(electoral_cycle_planned <= 1) %>% 
    mutate(failure = as.numeric(failure),
           regular = as.numeric(regular),
           opportunistic = as.numeric(opportunistic),
           minority_cabinet = as.numeric(minority_cabinet),
           largest_government_party = as.numeric(largest_government_party),
           single_party_cabinet = as.numeric(single_party_cabinet)) %>% 
    dplyr::select(poll_change,
                  poll_change_lag,
                  electoral_cycle_planned,
                  electoral_cycle_cabinet,
                  single_party_cabinet,
                  largest_government_party,
                  minority_cabinet,
                  gdp_change_lag,
                  gdp_change_lag_stand,
                  regular,
                  failure,
                  opportunistic,
                  diss_pm,
                  diss_gov,
                  vote_share_last) %>%
    pastecs::stat.desc(., basic = TRUE) %>%
    t() %>%
    as.data.frame() %>%
    tibble::rownames_to_column() %>%
    tbl_df %>%
    dplyr::select(rowname, nbr.val, mean, std.dev, min, median, max) %>%
    dplyr::rename(Variable = rowname, N = nbr.val, Mean = mean, SD = std.dev,
                  Min = min, Median = median, Max = max) %>%
    dplyr::mutate(Variable = nice.names)


print(desc_table)

# xtable::print.xtable(xtable::xtable(desc_table, 
#                                     digits=c(NA, NA, 0, rep(2, 5)),
#                                     caption="Descriptive statistics.",
#                                     label="tab:descriptives"),
#                      type="latex",
#                      #file="descriptives.tex",
#                      size = "footnotesize",
#                      html.table.attributes = "border=0",
#                      include.rownames=FALSE,
#                      sanitize.text.function=as.character,
#                      caption.placement="top")


## Table A2: Cases included in the analysis ----

summary_data <- dta_gov %>% 
    dplyr::filter(!is.na(gdp_change_lag),
                  !is.na(poll_change_lag)) %>%
    dplyr::group_by(country) %>%
    dplyr::summarise(`Electoral Cycles`=length(unique(election_id)),
                     `Mean Polls Per Cycle`=length(unique(poll_id)) / length(unique(election_id)),
                     `Date of First Poll`=format(min(polldate), "%d-%m-%Y"),
                     `Date of Last Poll`=format(max(polldate), "%d-%m-%Y")) %>%
    dplyr::rename(Country=country)

print(summary_data)

# xtable::print.xtable(xtable::xtable(summary_data, 
#                                     digits=c(NA, NA, 0, 1, NA, NA),
#                                     caption="Cases included in the analysis",
#                                     label="tab:summarydata",
#                                     align="llp{.15\\textwidth}p{.15\\textwidth}p{.15\\textwidth}p{.15\\textwidth}"),
#                      type="latex",
#                      #file="summarydata.tex",
#                      include.rownames=FALSE,
#                      sanitize.text.function=as.character,
#                      caption.placement="top")


### Regression models for models directly used in the paper

## Table A3 ----

screenreg(c(reg_base, reg_single_party_gov, reg_pm_diss),
       custom.coef.names = splitRecode(rownames(aggregate.matrix(
           lapply(c(reg_base, reg_single_party_gov, reg_pm_diss), extract))),
           by=":", recodes=nice.varnames),
       single.row = TRUE,
       include.aic=FALSE, include.bic=FALSE, include.variance=FALSE,
       custom.gof.names = c(NA, "N", "N (Party-Elections)", 
                            "N (Parties)", "N (Countries)"))
       
# texreg(c(reg_base, reg_single_party_gov, reg_pm_diss),
#        custom.coef.names = splitRecode(rownames(aggregate.matrix(
#            lapply(c(reg_base, reg_single_party_gov, reg_pm_diss), extract))),
#            by=":", recodes=nice.varnames),
#        single.row = TRUE,
#        include.aic=FALSE, include.bic=FALSE, include.variance=FALSE,
#        custom.gof.names = c(NA, "N", "N (Party-Elections)", "N (Parties)", "N (Countries)"),
#        #file="regressions_basic_paper_1-3.tex",
#        label = "table:m1-3",
#        float.pos = "h",
#        longtable=FALSE, use.packages=FALSE,
#        custom.model.names=paste("Model", 1:3),
#        custom.note = "%stars",
#        fontsize = "scriptsize",
#        caption.above= TRUE, caption= "Multilevel linear regression models of change in party support 
#        (coefficients with standard errors in parentheses).")

### Table A4 ----

screenreg(c(model_decades),
       custom.coef.names = splitRecode(rownames(aggregate.matrix(
           lapply(c(model_decades), extract))),
           by=":", recodes=nice.varnames),
       single.row = TRUE,
       no.space=TRUE,
       include.aic=FALSE, include.bic=FALSE, include.variance=FALSE,
       custom.gof.names = c(NA, "N", "N (Party-Elections)", 
                            "N (Parties)", "N (Countries)"))

# texreg(c(model_decades),
#        custom.coef.names = splitRecode(rownames(aggregate.matrix(
#            lapply(c(model_decades), extract))),
#            by=":", recodes=nice.varnames),
#        single.row = TRUE,
#        no.space=TRUE,
#        include.aic=FALSE, include.bic=FALSE, include.variance=FALSE,
#        custom.gof.names = c(NA, "N", "N (Party-Elections)", "N (Parties)", "N (Countries)"),
#        #file="regressions_basic_paper_4.tex",
#        label = "table:m4",
#        float.pos = "h",
#        longtable=FALSE, use.packages=FALSE,
#        custom.model.names=paste("Model", 4),
#        custom.note = "%stars",
#        fontsize = "scriptsize",
#        caption.above= TRUE, caption= "Multilevel linear regression models of change in party support 
#        (coefficients with standard errors in parentheses).")


## Note: Tables A5, A6, and A7 are part of the file models_plots_tables_nlme.R


## Table A8 ----

screenreg(c(reg_min_gov, reg_largest_gov_party),
       custom.coef.names = splitRecode(rownames(aggregate.matrix(
           lapply(c(reg_min_gov, reg_largest_gov_party), extract))),
           by=":", recodes=nice.varnames),
       single.row = F,
       include.aic = FALSE, include.bic = FALSE, include.variance = FALSE,
       custom.gof.names = c(NA, "N", "N (Party-Elections)",
                            "N (Parties)", "N (Countries)"))

# texreg(c(reg_min_gov, reg_largest_gov_party),
#        custom.coef.names = splitRecode(rownames(aggregate.matrix(
#            lapply(c(reg_min_gov, reg_largest_gov_party), extract))),
#            by=":", recodes=nice.varnames),
#        single.row = F,
#        include.aic = FALSE, include.bic = FALSE, include.variance = FALSE,
#        custom.gof.names = c(NA, "N", "N (Party-Elections)", "N (Parties)", "N (Countries)"),
#        #file="regressions_app_1-2.tex",
#        label = "table:app1-2",
#        float.pos = "h",
#        longtable = FALSE, use.packages = FALSE,
#        custom.model.names=paste0("Model A", 1:2),
#        custom.note = "%stars",
#        fontsize = "scriptsize",
#        caption.above = TRUE, caption = "Multilevel linear regression models of change in party support 
#        (coefficients with standard errors in parentheses).")


## Table A9 ----

screenreg(c(reg_lijphart, reg_gov_diss),
       custom.coef.names = splitRecode(rownames(aggregate.matrix(
           lapply(c(reg_lijphart, reg_gov_diss), extract))),
           by=":", recodes = nice.varnames),
       single.row = TRUE,
       include.aic = FALSE, include.bic = FALSE, include.variance = FALSE,
       custom.gof.names = c(NA, "N", "N (Party-Elections)", 
                            "N (Parties)", "N (Countries)"))

# texreg(c(reg_lijphart, reg_gov_diss),
#        custom.coef.names = splitRecode(rownames(aggregate.matrix(
#            lapply(c(reg_lijphart, reg_gov_diss), extract))),
#            by=":", recodes = nice.varnames),
#        single.row = TRUE,
#        include.aic = FALSE, include.bic = FALSE, include.variance = FALSE,
#        custom.gof.names = c(NA, "N", "N (Party-Elections)", "N (Parties)", "N (Countries)"),
#        #file="regressions_app_3-4.tex",
#        label = "table:app3-4",
#        float.pos = "h",
#        longtable = FALSE, use.packages = FALSE,
#        custom.model.names = paste0("Model A", 3:4),
#        custom.note = "%stars",
#        fontsize = "scriptsize",
#        caption.above = TRUE, caption = "Multilevel linear regression models of change in party support
#        (coefficients with standard errors in parentheses).")


## Table A10 ----

screenreg(c(reg_combined, reg_combined_subset),
       custom.coef.names = splitRecode(rownames(aggregate.matrix(
           lapply(c(reg_combined, reg_combined_subset), extract))),
           by=":", recodes = nice.varnames),
       single.row = FALSE,
       include.aic = FALSE, include.bic = FALSE, include.variance = FALSE,
       custom.gof.names = c(NA, "N", "N (Party-Elections)", 
                            "N (Parties)", "N (Countries)"))
       
# texreg(c(reg_combined, reg_combined_subset),
#        custom.coef.names = splitRecode(rownames(aggregate.matrix(
#            lapply(c(reg_combined, reg_combined_subset), extract))),
#            by=":", recodes = nice.varnames),
#        single.row = FALSE,
#        include.aic = FALSE, include.bic = FALSE, include.variance = FALSE,
#        custom.gof.names = c(NA, "N", "N (Party-Elections)", "N (Parties)", "N (Countries)"),
#        #file="regressions_app_5-6.tex",
#        label = "table:app5-6",
#        float.pos = "h",
#        longtable = FALSE, use.packages = FALSE,
#        custom.model.names=paste0("Model A", 5:6),
#        custom.note = "%stars",
#        fontsize = "scriptsize",
#        caption.above= TRUE, caption = "Multilevel linear regression models of change in party support 
#        (coefficients with standard errors in parentheses). Model 6 only includes countries with data for at least three cycles.")


## Table A11 ----

screenreg(c(reg_base, reg_base_lag, reg_base_percentage),
       custom.coef.names = splitRecode(rownames(aggregate.matrix(
           lapply(c(reg_base, reg_base_lag, reg_base_percentage), extract))),
           by=":", recodes=nice.varnames),
       single.row = TRUE,
       include.aic = FALSE, include.bic = FALSE, include.variance = FALSE,
       custom.gof.names = c(NA, "N", "N (Party-Elections)", 
                            "N (Parties)", "N (Countries)"))
       
# texreg(c(reg_base, reg_base_lag, reg_base_percentage),
#        custom.coef.names = splitRecode(rownames(aggregate.matrix(
#            lapply(c(reg_base, reg_base_lag, reg_base_percentage), extract))),
#            by=":", recodes=nice.varnames),
#        single.row = TRUE,
#        include.aic = FALSE, include.bic = FALSE, include.variance = FALSE,
#        custom.gof.names = c(NA, "N", "N (Party-Elections)", "N (Parties)", "N (Countries)"),
#        #file="regressions_app_7-9.tex",
#        label = "table:app7-9",
#        float.pos = "h",
#        longtable = FALSE, use.packages=FALSE,
#        custom.model.names = paste0("Model A", 7:9),
#        custom.note = "%stars",
#        fontsize = "scriptsize",
#        caption.above = TRUE, caption = "Multilevel linear regression models of change in party support 
#        (coefficients with standard errors in parentheses). 
#        Model A7 replicates Model 1, Model A8 adds the lagged DV, Model A9 uses the percentage change of a polling result compared to the previous election.")

## Table A12 ----

screenreg(c(reg_combined_gdp, reg_combined_gdp_stand),
       custom.coef.names = splitRecode(rownames(aggregate.matrix(
           lapply(c(reg_combined_gdp, reg_combined_gdp_stand), extract))),
           by=":", recodes=nice.varnames),
       single.row = TRUE,
       include.aic = FALSE, include.bic = FALSE, include.variance = FALSE,
       custom.gof.names = c(NA, "N", "N (Party-Elections)", "N (Parties)", "N (Countries)"))

# texreg(c(reg_combined_gdp, reg_combined_gdp_stand),
#        custom.coef.names = splitRecode(rownames(aggregate.matrix(
#            lapply(c(reg_combined_gdp, reg_combined_gdp_stand), extract))),
#            by=":", recodes=nice.varnames),
#        single.row = TRUE,
#        include.aic = FALSE, include.bic = FALSE, include.variance = FALSE,
#        custom.gof.names = c(NA, "N", "N (Party-Elections)", "N (Parties)", "N (Countries)"),
#        #file="regressions_app_10-11.tex",
#        label = "table:m1-3",
#        float.pos = "h",
#        longtable = FALSE, use.packages = FALSE,
#        custom.model.names = paste0("Model A", 10:11),
#        custom.note = "%stars",
#        fontsize = "scriptsize",
#        caption.above= TRUE, caption= "Multilevel linear regression models of change in party support 
#        (coefficients with standard errors in parentheses). 
#        Model A10 adds interaction of GDP Change and El. cycle to Model A5, Model A11 uses a measure of GDP change which is standardized by country and decade.")


## Table A13 ----

screenreg(c(reg_combined_elec_type, reg_combined, reg_combined_actual_cycle),
       custom.coef.names = splitRecode(rownames(aggregate.matrix(
           lapply(c(reg_combined_elec_type, reg_combined, reg_combined_actual_cycle), extract))),
           by=":", recodes=nice.varnames),
       single.row = TRUE,
       include.aic = FALSE, include.bic = FALSE, include.variance = FALSE,
       custom.gof.names = c(NA, "N", "N (Party-Elections)", 
                            "N (Parties)", "N (Countries)"))
       
# texreg(c(reg_combined_elec_type, reg_combined, reg_combined_actual_cycle),
#        custom.coef.names = splitRecode(rownames(aggregate.matrix(
#            lapply(c(reg_combined_elec_type, reg_combined, reg_combined_actual_cycle), extract))),
#            by=":", recodes=nice.varnames),
#        single.row = TRUE,
#        include.aic = FALSE, include.bic = FALSE, include.variance = FALSE,
#        custom.gof.names = c(NA, "N", "N (Party-Elections)", "N (Parties)", "N (Countries)"),
#        #file="regressions_app_12-14.tex",
#        label = "table:m11-13",
#        float.pos = "h",
#        longtable=FALSE, use.packages=FALSE,
#        custom.model.names = paste0("Model A", 12:14),
#        custom.note = "%stars",
#        fontsize = "scriptsize",
#        caption.above= TRUE, caption= "Multilevel linear regression models of change in party support 
#        (coefficients with standard errors in parentheses). 
#        Model A12 includes the type of election, 
#        Model A13 reproduces Model A5 to compare the coefficients to Model A14 that uses the planned length of the electoral cycle for each country.")


